function [Emax,Emin,Phi,eEmax,eEmin,n0,flc,fls,u0,v0,eu0,ev0,...
    omeg,eomeg] = gs_strain(em,nm,e,n,ve,vn,eve,evn,npe,d0,Options)
%gs_strain - Strain calculation on a single grid point
%
%   [Emax,Emin,Phi,eEmax,eEmin,n0,flc,fls,u0,v0,eu0,ev0,...
%       omeg,eomeg] = gs_strain(em,nm,e,n,ve,vn,eve,evn,npe,d0,Options)
%
% This function calculates the strain rate in the grid node (em,nm), where
% 'e' indicates 'East direction' (x-axis) and 'n' is 'North direction'
% (y-axis), on the basis of velocities (or displacements) (ve,vn) at the
% measurement points whose coordinates are (e,n).
% If npe(k) is true, the point (e(k),n(k)) is effectively considered in
% the calculations. If npe(k) is false, the corresponding experimental
% point is not considered.
% The inputs eve, evn are the errors on ve, vn respectively.
%
% Output data are:
%
%   Emax   : Maximum eigenvalue (max strain);
%   Emin   : Minimum eigenvalue (min strain);
%   Phi    : Azimuth of direction of maximum strain;
%   eEmax  : Error on Emax; 
%   eEmin  : Error on Emin;
%   n0     : # of experimental points whose distance from (em,nm) is lower
%            than d0 ('critical distance')
%   flc    : Flag on calculations. It is 1 if the calculation has been made, 
%            0 otherwise. If flc is 0, all the others output data are NaN
%   fls    : Flag on significance. It is 2 if all the three sectors with
%            angular width 120� around the point (em,nm) contain at least 
%            one experimental point whose distance by (em,nm) is lower than
%            d0 (HIGH SIGNIFICANCE GRID POINT CASE).
%            It is 1 if two sectors contain at least one point as above
%            (MID SIGNIFICANCE GRID POINT CASE).
%            It is 0 othewise (LOW SIGNIFICANCE GRID POINT CASE).
%            Obvioulsly, fls is always 0 if flc = 0;
%   u0,v0  : Translation along East and North axis respectively, obtained
%            as 'complementary data' of the least square computation;
%   eu0,ev0: corresponding errors;
%   omeg   : rotation;
%   eomeg  : corresponding error.
%
% This function is called by GridStrain (strain calculation on a grid), 
% but it can also be used as independent function; the use of a grid is 
% not compulsory. 
%
%   [Emax,Emin,Phi,eEmax,eEmin,n0,flc,fls,u0,v0,eu0,ev0,...
%       omeg,eomeg] = gs_strain(em,nm,e,n,ve,vn,eve,evn,npe,d0)
%
% See also gs_itera.

% A. Pesci &  G. Teza, 2004, 2005, 2008, 2011, 2022.  

format long e; 
warning off all;

if strcmpi(Options.Function,'exp')
    optfun = 1;
elseif strcmpi(Options.Function,'gaussian')
    optfun = 2;
else
    optfun = 3;
end

le = length(e);
N  = le*2;

OBS = zeros(N,1);  % Observation vector initializing
OBS(1:2:N) = ve;
OBS(2:2:N) = vn;

de = e-em;
dn = n-nm;
vd = [0; 1]; % north direction unit vector

% Design matrix generation; it is a N-by-6 matrix, where
% N = 2*M (M is the number of observed points). This matrix 
% is related to the point (em,nm).
%
%   |1 0 de1 dn1  0   0 |         |u1|
%   |0 1  0   0  de1 dn1|         |v1|
%   |1 0 de2 dn2  0   0 | |u0 |   |u2|
%   |0 1  0   0  de2 dn2| |v0 |   |v2|
%   |. .  .   .   .   . | |L11|   |..| 
%   |. .  .   .   .   . | |L12| = |..|
%   |. .  .   .   .   . | |L21|   |..| 
%   |. .  .   .   .   . | |L22|   |..|
%   |1 0 deM dnM  0   0 |         |uM| 
%   |0 1  0   0  deM dnM|         |vM|
  
A = zeros(N,6);
A(1:2:N-1,1) = 1;
A(2:2:N,2)   = 1;
A(1:2:N-1,3) = de;
A(1:2:N-1,4) = dn;
A(2:2:N,5)   = de;
A(2:2:N,6)   = dn;

% Generation of error's vector...
EV = zeros(N,1);
EV(1:2:N) = eve;
EV(2:2:N) = evn;
EV = EV.^2;

%... and its weighting using an exponential
def = d0;         % definition of "significance distance" (for a flag)
dis = sqrt(de.^2+dn.^2);        % distances vector
n0  = sum(dis <= d0);  % computation of experimental points 
                  % whose distance from (em,nm) is lower than d0.  
% error scaler:
if optfun == 1 
    scalp = exp(2*dis/d0); 
elseif optfun == 2  
    scalp = exp((dis/d0).^2);
else
    scalp = 1./(1+(dis/d0).^2);
end
scal  = zeros(N,1);
scal(1:2:N-1) = scalp;
scal(2:2:N)   = scalp;
EVS   = EV.*scal; % In this way, the errors on each experimental point 
                  % are distance-weighted (if a point is far from (em,nm), 
                  % its is magnified)
                       
V = diag(EVS);    % Standard deviations

% implementation of least squares method
flc = 0;
fls = 0;
if nonaninf(V)
    
    P = V\eye(size(V));         % accurate computation of P = inv(V)
    B = A'*P*A;
    
    if nonaninf(B)
        
        Binv = B\eye(size(B));  % accurate computation of Binv = inv(B)
        GRAD = Binv*A'*P*OBS;   % Estimated solution (U, V, gradients)
                % GRAD elements are u0, v0, L11, L12, L21 and L22
        OBS1 = A*GRAD;      % Expected observation from estimated solution
        OBR  = OBS1-OBS;    % Observation residual
        SIG2 = (OBR'*P*OBR)/(N-6);   % 6 elements are estimated
        DISP = Binv*SIG2;   % Dispersion of solution, i.e. variance
        eDISP = diag(DISP);
        ST = [GRAD(3), (GRAD(4)+GRAD(5))/2; (GRAD(4)+GRAD(5))/2, GRAD(6)];
                            % strain tensor
        eST = zeros(2);
        eST(1,1) = sqrt(eDISP(3));
        eST(1,2) = sqrt(eDISP(4)+eDISP(5))/2;
        eST(2,1) = eST(1,2);
        eST(2,2) = sqrt(eDISP(6));
        
        u0 = GRAD(1); 
        v0 = GRAD(2);
        eu0 = sqrt(eDISP(1)); 
        ev0 = sqrt(eDISP(2)); 
        
        omeg = (GRAD(4)-GRAD(5))/2;   % rotation
        eomeg = eST(1,2);             % corresponding uncertainty 
        
        try %#ok<TRYNC>
            
            [TM,EVM] = eig(ST);  % EVM is a diagonal matrix, whose elements
                                 % are the eigenvalues; TM is the base transformation 
                                 % matrix. It is inv(TM)*ST*TM = EVM.
            Emax = EVM(2,2);     % Eig provides eigenvnalues sorted in 
            Emin = EVM(1,1);     % ascending order! (the maximum is the 2nd!)

            % Computation of Phi starting from normalized eigenvectors 
            Phi  = acos(dot([0;1],TM(:,2)));  % OK if east > 0, north > 0
            if (TM(1,2) < 0)&&(TM(2,2) > 0)       
                Phi = -Phi;                   % OK if east < 0, north > 0
            elseif (TM(1,2) < 0)&&(TM(2,2) < 0)
                Phi = pi-Phi;                 % OK if east < 0, north < 0
            elseif (TM(1,2) > 0)&&(TM(2,2) < 0)
                Phi = Phi-pi;                 % OK if east > 0, north < 0
            end
            
            eEVM  = TM\(eST*TM);  % Error matrix transformation (inv(TM)*eST*TM)
            eEmax = eEVM(2,2);
            eEmin = eEVM(1,1);
%             if abs(eEVM(1,2)) > 0.5*max(eEVM(1,1),eEVM(2,2))
%                 eEmax = 3*eEmax;    % Error maximization if the covariance
%                 eEmin = 3*eEmax;    % is a significant part of STD
%             end
            
            flc = 1;  % computation ok
           
            %--- analysis of calculation significance
            if ~Options.HSFourQuadrants
                ntv = zeros(1,3);
            else
                ntv = zeros(1,4);
            end
            for kk = 1:le
                if npe(kk) == 1 
                    vg  = [de(kk); dn(kk)];
                    dvg = dis(kk);    % already defined, see above
                    if dvg <= def
                        if dvg > 100*eps
                            vgn = vg/dvg;
                            aze = acos(dot(vd,vgn));
                            if de(kk) < 0
                                aze = 2*pi-aze; 
                            end
                            if ~Options.HSFourQuadrants
                                if aze < 2*pi/3     % first  120� sector
                                    ntv(1) = ntv(1)+1;
                                elseif aze < 4*pi/3 % second 120� sector
                                    ntv(2) = ntv(2)+1;
                                else                % third  120� sector
                                    ntv(3) = ntv(3)+1;
                                end
                            else
                                if aze < pi/2
                                    ntv(1) = ntv(1)+1;
                                elseif aze < pi
                                    ntv(2) = ntv(2)+1;
                                elseif aze < 3*pi/2
                                    ntv(3) = ntv(3)+1;
                                else
                                    ntv(4) = ntv(4)+1;
                                end
                            end
                        end
                    end
                end
            end
            if sum(logical(ntv)) == numel(ntv)
                fls = 2;
            elseif sum(logical(ntv)) == numel(ntv)-1 
                fls = 1;
            end
            
        end
        
    end
    
else
    
    u0 = NaN; v0 = NaN; eu0 = NaN; ev0 = NaN;
    
end

if flc == 0
    Emax  = NaN; Emin  = NaN;
    Phi   = NaN;
    eEmax = NaN; eEmin = NaN;
end